{
 "metadata": {
  "name": "",
  "signature": "sha256:4e7d7d565dabddb1455a5227da85112f04a86bd1f76179c5f08615a26caaa5fa"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Zhenfeng Liang"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# MTH9879 Homework 10\n",
      "\n",
      "Assigned: April 28, 2015.\n",
      "Due: May 12, 2015 by 6pm. \n",
      "\n",
      "Late homework **will not be accepted**.\n",
      "\n",
      "$$\n",
      "\\newcommand{\\supp}{\\mathrm{supp}}\n",
      "\\newcommand{\\E}{\\mathbb{E}}\n",
      "\\newcommand{\\Eof}[1]{\\mathbb{E}\\left[ #1 \\right]}\n",
      "\\def\\Cov{{ \\mbox{Cov} }}\n",
      "\\def\\Var{{ \\mbox{Var} }}\n",
      "\\newcommand{\\1}{\\mathbf{1} }\n",
      "\\newcommand{\\PP}{\\mathbb{P} }\n",
      "%\\newcommand{\\Pr}{\\mathrm{Pr} }\n",
      "\\newcommand{\\QQ}{\\mathbb{Q} }\n",
      "\\newcommand{\\RR}{\\mathbb{R} }\n",
      "\\newcommand{\\DD}{\\mathbb{D} }\n",
      "\\newcommand{\\HH}{\\mathbb{H} }\n",
      "\\newcommand{\\spn}{\\mathrm{span} }\n",
      "\\newcommand{\\cov}{\\mathrm{cov} }\n",
      "\\newcommand{\\sgn}{\\mathrm{sgn} }\n",
      "\\newcommand{\\HS}{\\mathcal{L}_{\\mathrm{HS}} }\n",
      "%\\newcommand{\\HS}{\\mathrm{HS} }\n",
      "\\newcommand{\\trace}{\\mathrm{trace} }\n",
      "\\newcommand{\\LL}{\\mathcal{L} }\n",
      "%\\newcommand{\\LL}{\\mathrm{L} }\n",
      "\\newcommand{\\s}{\\mathcal{S} }\n",
      "\\newcommand{\\ee}{\\mathcal{E} }\n",
      "\\newcommand{\\ff}{\\mathcal{F} }\n",
      "\\newcommand{\\hh}{\\mathcal{H} }\n",
      "\\newcommand{\\bb}{\\mathcal{B} }\n",
      "\\newcommand{\\dd}{\\mathcal{D} }\n",
      "\\newcommand{\\g}{\\mathcal{G} }\n",
      "\\newcommand{\\p}{\\partial}\n",
      "\\newcommand{\\half}{\\frac{1}{2} }\n",
      "\\newcommand{\\T}{\\mathcal{T} }\n",
      "\\newcommand{\\bi}{\\begin{itemize}}\n",
      "\\newcommand{\\ei}{\\end{itemize}}\n",
      "\\newcommand{\\beq}{\\begin{equation}}\n",
      "\\newcommand{\\eeq}{\\end{equation}}\n",
      "\\newcommand{\\beas}{\\begin{eqnarray*}}\n",
      "\\newcommand{\\eeas}{\\end{eqnarray*}}\n",
      "\\newcommand{\\cO}{\\mathcal{O}}\n",
      "\\newcommand{\\cC}{\\mathcal{C}}\n",
      "\\newcommand{\\cF}{\\mathcal{F}}\n",
      "\\newcommand{\\cL}{\\mathcal{L}}\n",
      "\\newcommand{\\BS}{\\text{BS}}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<font color = \"red\">Homework is to be done by each student individually.  To receive full credit, you must email a completed copy of this iPython notebook to Yu Gan (yugan323@gmail.com), Fubo Shi (fubo.shi.baruch@gmail.com), and Tai-Ho Wang (tai-ho.wang@baruch.cuny.edu) by the due date and time.  All R-code must run correctly and solutions must be written up neatly in Markdown/LaTeX format.\n",
      "\n",
      "<font color=\"blue\">If you encounter problems with Markdown/LaTeX or iPython notebook, please contact your TAs Yu Gan and/or Fubo Shi.\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Set up R environment"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext rmagic"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Numerical solution of Fredholm equation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 1. (12 points)\n",
      "Consider the Fredholm equation (7) in Lecture 10:\n",
      "\n",
      "$$\n",
      "\\int_0^T{G}(|t-s|)\\,dx_s =\\text{const.}\\qquad\\text{for all $t\\in [0,T]$.}\n",
      "$$\n",
      "\n",
      "Define $t_i:={i}\\frac T{N}$. One way to discretize this integral equation is as follows:\n",
      "\n",
      "$$\n",
      "\\sum_{j=1}^N G_{ij}\\,v_j=\\text{const.}\n",
      "$$\n",
      "\n",
      "where $v_j=x_{t_{j}}-x_{t_{j-1}}$ and\n",
      "\n",
      "<a name=\"eq:1\"></a>(1)\n",
      "$$\n",
      "G_{ij}=\n",
      "\\int_{t_{i-1}}^{t_i}\\,\\int_{t_{j-1}}^{t_j}\\,G(|t-s|)\\,ds\\,dt.\n",
      "$$\n",
      "\n",
      "(a) Derive an explicit expression for the $G_{ij}$ in the case\n",
      "$$\n",
      "G(\\tau)=\\frac {1}{(1+\\tau)^2}.\n",
      "$$\n",
      "You will need to compute separate expressions for the cases $i \\neq j$ and $i=j$.\n",
      "\n",
      "(b) With $N=512$, use matrix algebra to solve equation [(1)](#eq:1) numerically with this decay kernel.\n",
      "\n",
      "(c) Plot your results - you should end up with something like the top figure in Figure 2 of Lecture 10.\n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Solution\n",
      "\n",
      "### (a)\n",
      "\n",
      "We need to consider three cases, $i = j$, $i < j$ and $i > j$.\n",
      "\n",
      "Case 1, $i < j$\n",
      "\n",
      "\\begin{eqnarray*}\n",
      "G_{ij}&=&\n",
      "\\int_{t_{i-1}}^{t_i}\\,\\int_{t_{j-1}}^{t_j}\\,G(|t-s|)\\,ds\\,dt. \\\\\n",
      "&=& \\int_{t_{i-1}}^{t_i}\\,\\int_{t_{j-1}}^{t_j}\\,\\frac{1}{(1+s-t)^2}\\,ds\\,dt.\n",
      "\\end{eqnarray*}\n",
      "\n",
      "Case 2, $i > j$\n",
      "\n",
      "\\begin{eqnarray*}\n",
      "G_{ij}&=&\n",
      "\\int_{t_{i-1}}^{t_i}\\,\\int_{t_{j-1}}^{t_j}\\,G(|t-s|)\\,ds\\,dt. \\\\\n",
      "&=& \\int_{t_{i-1}}^{t_i}\\,\\int_{t_{j-1}}^{t_j}\\,\\frac{1}{(1+t-s)^2}\\,ds\\,dt.\n",
      "\\end{eqnarray*}\n",
      "\n",
      "Case 3, $i = j$\n",
      "\n",
      "\\begin{eqnarray*}\n",
      "G_{ij}&=&\n",
      "\\int_{t_{i-1}}^{t_i}\\,\\int_{t_{j-1}}^{t_j}\\,G(|t-s|)\\,ds\\,dt. \\\\\n",
      "&=& \\int_{t_{i-1}}^{t_i}\\,\\int_{t_{j-1}}^{t}\\,\\frac{1}{(1+t-s)^2}\\,ds + \\,\\int_{t}^{t_j}\\,\\frac{1}{(1+s-t)^2}\\,ds\\,dt.\n",
      "\\end{eqnarray*}\n",
      "\n",
      "With some calculus, we got,\n",
      "\n",
      "$$\n",
      "G_{ij} = \\begin{cases} \\ln(1 + t_j - t_i) - \\ln(1 + t_j - t_{i-1}) + \\ln(1 + t_{j-1} - t_{i - 1}) - \\ln(1 + t_{j-1} - t_{i}), & \\mbox{if } i < j \\\\ \n",
      "\\ln(1 + t_i - t_j) - \\ln(1 + t_{i-1} - t_{j}) + \\ln(1 + t_{i-1} - t_{j - 1}) - \\ln(1 + t_{i} - t_{j-1}), & \\mbox{if } i > j \\\\ \n",
      "2\\,(t_i - t_{i-1}) -2\\,\\ln(1 + t_i - t_{i-1})  & \\mbox{if } j = i\\end{cases}\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "\n",
      "N = 512\n",
      "T = 1\n",
      "\n",
      "cal_G <- function(i, j){\n",
      "\n",
      "    t_i = i * T / N\n",
      "    t_j = j * T / N\n",
      "    t_i_minus_1 = (i - 1) * T / N\n",
      "    t_j_minus_1 = (j - 1) * T / N\n",
      "    \n",
      "    if(i < j)\n",
      "    {\n",
      "        result <- log(1 + t_j - t_i) - log(1 + t_j - t_i_minus_1) + log(1 + t_j_minus_1 - t_i_minus_1) - log(1 + t_j_minus_1 - t_i)\n",
      "    }\n",
      "    else if(i > j)\n",
      "    {\n",
      "        result <- log(1 + t_i - t_j) - log(1 + t_i_minus_1 - t_j) + log(1 + t_i_minus_1 - t_j_minus_1) - log(1 - t_j_minus_1 + t_i)\n",
      "    \n",
      "    }\n",
      "    else\n",
      "    {\n",
      "        result <- 2 * (t_i - t_i_minus_1) - 2 * log(1 + t_i - t_i_minus_1)\n",
      "    }\n",
      "        \n",
      "    return (result)\n",
      "} \n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "\n",
      "G <- matrix(, nrow = N, ncol = N)\n",
      "\n",
      "for (i in 1 : N)\n",
      "{\n",
      "    for (j in 1 : N)\n",
      "    {\n",
      "        G[i, j] <- cal_G(i,j)\n",
      "    }    \n",
      "}\n",
      "\n",
      "const <- rep(1/100, N)\n",
      "\n",
      "# Solve for v\n",
      "v <- solve(G, const)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "\n",
      "## (c)\n",
      "\n",
      "# Construct time\n",
      "t <- seq(T/N,T,by=T/N)\n",
      "\n",
      "# plot\n",
      "plot(t, v, col='red', type='l', main='v and time', ylab='v', xlab='time')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAIAAADytinCAAAeUUlEQVR4nO3de1BU9/3/8c9hF1Dc\nVVGMggVpqKgjaqvYqlGDyqjRxgka4qXTVjtiWnU6NZpM0+lE63QmTaPVZDodEzVTO3ZixwYN1rs1\n5NKBZJJ+46WF1luCRlBEFNQsAXZ/f6w/K6gre+R8zvscno+/lt2Fz/vA+Pq8PHsWjFAopAAA8sTY\nPQAA4O4IaAAQioAGAKEIaAAQioAGAKEIaAAQioAGAKEIaAAQioAGAKEIaAAQioAGAKEIaAAQioAG\nAKEIaAAQioAGAKEIaAAQioAGAKEIaAAQioAGAKEIaAAQioAGAKEIaAAQioAGAKEIaHQ4x48fNwxj\n4MCB932mYRiGYTQ1NWmYCrgTAQ38T3Z2tmEYpaWldg8CKKWU1+4BALnq6+uVUh6Px+5B0EHRoCHI\nD3/4Q8MwXnzxxfCHzz77rGEYy5cvb/W0Q4cOjR492u/3JyYmTpo06ZNPPlFKVVVVGYaRlJS0Y8eO\nrKwsv9//+OOPX7hwIfwpVVVVeXl5iYmJw4YNu1dBzs7ODn+p0aNHb9u2LRAI+P3+7t27G4Zx8uRJ\nwzAyMjLWrVv3ta99rV+/fq+88kppaenw4cN9Pl9OTs7p06fDX+T06dMzZsxISkrq27fvD37wg+rq\naiu+UegoQoAYu3fvVkqNGzcu/OGwYcOUUiUlJbc/p6KionPnzh6PJycnZ+TIkUqp1NTUYDBYWVmp\nlPJ4PH6/f/jw4TExMUqphQsXhkKhxsbGwYMHh5/5rW99K/zQgAED7lw9PT1dKbVy5cozZ858+eWX\n4S8YCoVOnDgR/vcSFxc3YMCAW7czMjK6dOmilMrLywuFQnV1db179/Z4PN/97nfHjh2rlMrKymps\nbNTxvYMbEdAQpKGhITEx0ePxXL58OVx+w+F7+3OKi4unTJmyevXq8PM7deqklLp48WI4oJVSn376\naSgU2rJli1JqyJAhoVDorbfeUkoNHjz4+vXrwWBw0aJFdw3oUCg0YsSIW1vCXQO6rKwsFAqNHz9e\nKTV37txgMHj48GGlVEZGRigUWr9+vVLqRz/6UXV1dXV1dXZ2tlJqx44dln/j4FKc4oAgcXFxM2fO\nbG5uPnjw4KFDh5RS+fn5hmHc/pxHH3301Vdf9Xg8eXl5/fv3DwQCSqnm5ubwo+GTGEqpcLm+ceOG\nUurIkSNKqVmzZiUkJBiGMW/ePHPj9e3bN3ztR2pqqlJq4sSJhmGEb4ev9Dh+/LhS6o033ujVq1ev\nXr0+/vhjpdS//vUvc8sBvEgIWWbPnr158+a9e/eGP3zqqadaPaGkpOTRRx+NjY3Ny8v7xS9+8fOf\n//zKlSu3Hg2fvlBK3R7rwWDw9ntMv+jn9XojfKiUamhoUEqtWLFiypQpt+5MS0sztxxAg4YsEyZM\n6NWr1969ew8ePJiWlvbtb3+71RMKCwsbGxuXLl26devWyZMn357O95KVlRX+xPBZizfffDPy82/1\n8WgNGjRIKRUIBHJzc3Nzczt16nTp0iUuAoFpNGjI4vV6Z82atWHDBqXU8uXLW53fUEr17t1bKfXa\na6+VlZV9+OGHhmGEQqHm5uYIOZiXl/fwww8fO3Zs0KBBPXv2/Oc//3mvZ4Zf8fvd737XvXv3jIyM\naIdfuHDhSy+99Ic//OHy5cvBYHD79u1+v//f//53tF8HCKNBQ5zZs2eHb+Tn59/56E9+8pOZM2c2\nNjaePHly/fr1jz76qFIqfPnHvcTFxb3zzjuPPfZYbW1tIBBYu3btvZ75s5/9rHfv3oWFhceOHTMx\nea9evT744INJkybt3r17//7906dP/+CDD5KTk018KUApZYRCIbtnAADcBQ0aAIQioAFAKAIaAIQi\noAFAKAIaAIQioAFAKAIaAIQioAFAKAIaAIQioAFAKAIaAIQioAFAKAIaAIQioAFAKAIaAIQioAFA\nKAIaAIQioAFAKAIaAIQioAFAKAIaAIQioAFAKK/dA0SnpqamsLAwFArZPQgAKKVUfHz8vHnzYmNj\nrfjiDmvQf//734uLi+2eAgBu2rhxY0VFhUVf3GENWin1yCOPLFq0yO4pAEAppT766CPrvrjDGjQA\ndBwENAAIRUADgFAENAAIRUADgFAENAAIRUADgFAENAAIRUADgFAENAAI5eqAPnDA7gkAOFBjo3rn\nHbuHUMrlAf3ii3ZPAMCBLl9Wr71m9xBKuTygAcDJCGgAaCkUUoZh9xBKEdAAIBYBDQAt0aABAJER\n0ADQEg0aAIQioAEAkRHQANASDRoAEBkBDQAt0aABAJER0ADQEg0aABAZAQ0ALdGgAUAoAhoAEBkB\nDQAt0aABAJER0ADQEg0aABAZAQ0ALdGgAUAoAhoAEBkBDQAt0aABAJER0ADQEg0aABAZAQ0ALdGg\nAQCREdAA0BINGgCEIqABAJER0ADQEg0aABAZAQ0ALdGgAQCREdAA0BINGgCEIqABAJER0ADQEg0a\nABAZAQ0ALdGgAQCREdAA0BINGgAQGQENAC3RoAFAKAIaABAZAQ0ALdGgAQCREdAA0BINGgAQmasD\n2jBUKGT3EACchgatAwENwISOGdDBYLCuri4YDOpcFAAcSkdABwKBlStXZmZmxsfHd+vWLS4urn//\n/qtWrWpoaLB2YRo0ABM6VINetGhRSUnJxo0bq6qqvvrqq4sXL27ZsuXo0aOLFy/WsDoAOJRXwxpF\nRUVlZWXJycnhD3v06DFmzJitW7f269fP2oVp0ABM6FANOj09fd++fa3u3L9/f2pqqobVAcChdDTo\nTZs2zZgxY82aNVlZWX6/v76+vqysrKampqioyNqFadAATBDToHUEdHZ2dkVFRXFx8enTp2traxMT\nEwsKCnJycrxeHasDgENpikiv15ubmxsMBq9du+bz+WJitFzeR4MGYIKYBs1ldgDQUocKaC6zAwAT\nuMwOAFrqUA2ay+wAwAQuswOAlsQ0aC6zAwCh5F5mt3379tdff73VnSdPnhwwYEBbX12kQQMwoUM1\n6EAg8OKLL7755ptnzpxpamryeDxf//rXv/e97z3//PPx8fH3+qz8/Pz8/PxWdy5btqyysrKtCxPQ\nAEwQE9BcZgcAQnGZHQC0JCY3uMwOAO4g4xQHl9kBQEticoPL7ADgDh2nQav/f5mdnrX+hwYNwAQx\nuaH1r3oDgDN0nAZdXl5+r4cGDhxo4cI0aAAmiLkOWkdAP/PMM3v37k1ISEhMTGz10Llz5yxcmIAG\nYIKY3NAR0Hv27CkoKIiPj//973+vYTkAeFAyGrSmc9Bz5sxJT0/Xs9b/0KABmCAmNzRdxTFp0qRJ\nkybpWQsAHlSHatD2oEEDMEFMbrg6oAHAHBq05WjQAEwQc5kdAQ0ALYnJDVcHNACYQ4O2HA0agAli\ncsPVAQ0A5tCgLUeDBmCCmNxwdUADgDk0aMvRoAGYICY3XB3QAGAODdpyNGgAJvBGFR0IaAAmiMkN\nVwc0AJhDg7YcDRqACWJyw9UBDQDm0KAtR4MGYIKY3HB1QAOAOTRoy9GgAZjAZXY6ENAATBCTG64O\naAAwhwZtORo0ABPE5IarAxoAzKFBW44GDcAEMbnh6oAGAHNo0JajQQMwQUxuuDqgAcAcGrTlaNAA\nTOCNKjoQ0ABMEJMbrg5oADCHBm05GjQAE8TkhqsDGgDMoUFbjgYNwAQxueHqgAYAc2jQlqNBAzCB\ny+x0IKABmCAmN1wd0ABgDg3acjRoACaIyQ1XBzQAmEODthwNGoAJYnLD1QENAObQoC1HgwZggpjc\ncHVAA4A5NGjL0aABmMAbVXQgoAGYICY3XB3QAGAODdpyNGgAJojJDVcHNACYQ4O2HA0agAlicsPV\nAQ0A5tCgLUeDBmACl9npQEADMEFMbrg6oAHAHBq05WjQAEwQkxuuDmgAMIcGbTkaNAATxOSGqwMa\nAMyhQVuOBg3ABDG54eqABgBzaNCWo0EDMIE3quhAQAMwQUxuuDqgAcAcGrTlaNAATBCTG64OaAAw\nhwZtORo0ABPE5IarAxoAzKFBW44GDcAELrPTgYAGYIKY3NAa0MFgsK6uLhgM6lwUAKLWcRp0IBBY\nuXJlZmZmfHx8t27d4uLi+vfvv2rVqoaGBmsXpkEDMEFMbugI6EWLFpWUlGzcuLGqquqrr766ePHi\nli1bjh49unjxYg2rA0DUZDRor4Y1ioqKysrKkpOTwx/26NFjzJgxW7du7devn7UL06ABmCAmN3Q0\n6PT09H379rW6c//+/ampqRpWB4CodZwGvWnTphkzZqxZsyYrK8vv99fX15eVldXU1BQVFVm7MA0a\ngAlicuMuAT1nzpy5c+dOmTKlU6dO7bJGdnZ2RUVFcXHx6dOna2trExMTCwoKcnJyvF4d2wMARE1s\ngx4xYsRvf/vb+fPnP/HEE3PmzJk4cWJsbOyDLuP15ubmBoPBa9eu+Xy+mBgtl/fRoAGYIPmNKs8+\n++w//vGPsrKy0aNHr1+/Pj09/cc//vGDrMFldgCcRExu3LPJdu/ePS0t7Rvf+IZhGO+///6DrMFl\ndgAcRkaDvsspjtdff3337t2HDx/OysrKy8s7fPhwZmbmg6xh7jK7v/71r7/5zW9a3Xnu3LkhQ4a0\ndWEaNAATxOTGXQK6sLDwiSee2LBhw61IfUDhy+wWLFhw+533vczuySeffPLJJ1vduWzZssrKynaZ\nCgDuSWyDvvOa5QfEZXYAnERMbui40I3L7AA4jNgGbckyXm9ubm74dnV1dWxsrI50pkEDMEHyZXbt\n7vHHHz979qxS6ty5c6NHj+7Tp89DDz00ceLE8+fPW742AQ0gWmJyQ0dAHzx48Pr160qp5cuXDxw4\nsL6+/vr16yNHjly6dKm1C8vYAwE4j4z00HoW+OOPP96zZ09CQoJS6vnnn7f8t9kpQTshAMcQkxua\n/qJKZWVlU1PT4MGDT506Fb7n2LFjPp/P2lVl7IEAnEdGeuho0OPHj58/f/6FCxc6d+584sSJadOm\nFRcX5+XlvfDCC5avLWYnBOAYYnJDR0AfOHBAKdXY2FhRUVFVVaWU6ty5844dO3JycqxdWMYeCMB5\nZKSHvnPQsbGxGRkZGRkZSqnvfOc7mlYVsxMCcAwxuaH1r3rrJmMPBOA8MtLD7QEtZicE4Bgd6o0q\ndiKgAURLTG64OqBl7IEAnEdGerg6oJWgnRCAY4jJDVcHtIw9EIDzyEgPVwe0ErQTAnAMMbnh6oCW\nsQcCcB4Z6eH2gBazEwJwDC6z04SABhAtMbnh6oCWsQcCcB4Z6eHqgFaCdkIAjiEmN1wd0DL2QADO\nIyM9XB3QStBOCMAxxOSGqwNaxh4IwHlkpIerA1oJ2gkBIFquDmgZeyAAh+E6aB14owoAEwhoTQho\nAI7l6oCWsQcCcBgatCY0aACO5eqAlrEHAnAYGrQmNGgAjuXqgJaxBwJwGBq0DlxmB8AEAloTAhqA\nY7k6oGXsgQAchgatCQ0agGO5OqBl7IEAHIYGrQkNGoBjuTqgZeyBAByGBq0JDRqAY7k6oGXsgQAc\nhgatA29UAWACAa0JAQ3AsVwd0DL2QAAOQ4PWhAYNwLFcHdAy9kAADkOD1oQGDcCxXB3QMvZAAA5D\ng9aBy+wAmEBAa0JAA3AsVwe0jD0QgMPQoDWhQQNwLFcHtIw9EIDD0KA1oUEDcCxXB7SMPRCAw9Cg\nNaFBA3AsVwe0jD0QgMPQoHXgjSoATCCgNSGgATiWqwNaxh4IwGFo0JrQoAE4lqsDWsYeCMBhaNCa\n0KABOJarA1rGHgjAYWjQOnCZHQATCGhNCGgAjuXqgJaxBwJwGBq0JjRoAI7l6oCWsQcCcBgatCZv\nvKHefdfuIQA4x/r16v/+z+4hbnJ1QBuGqq5WFy/aPQcA56isVLW1NGgtgkEVCNg9BADnCARUMGj3\nEDe5OqANQzU3E9AAohAIqOZmGrT1CGgA0SKg9eEUB4CoNDRwikOLcINuaLB7DgDOQYPWJxgkoAFE\ngRcJNeEcNIBo0aD1IaABRCUc0DK4OqBp0ACiRYPWhIAGEK2GBgJaFy6zAxAVXiTUhAYNIFqc4tCH\ny+wARIUGrUl4D6RBA2i7hoYO+vugg8FgXV1dUPPuREADaKNQSFRi6AjoQCCwcuXKzMzM+Pj4bt26\nxcXF9e/ff9WqVQ1Wn3wI74E1NSonx9qFALhAXZ0aN+7m38mT0aC9GtZYtGhRVVXVxo0bs7Kyunbt\nWl9fX15evmbNmsWLF2/evNny5T//XH32mZoxQ8XGquvXVXOz8vlUbKzl6yqlYmOVz6djIW28XuX3\n2z0E7lBba/cE7e3qVU0ngr/8UgUCKi5Oxcaqq1fVRx/pWLTNdAR0UVFRWVlZcnJy+MMePXqMGTNm\n69at/fr1i/xZW7dubXXnkSNH0tLS2rpweA9sblZLlqgFC1RGRvv/DdmmJlVf385f0wqNjeraNbuH\ngDW8XtXUdP87W91z1yf4fCKa463Zwjfu/LC9WkL4C3o8KhhUHo8qKVGvvKL27FGqIzXo9PT0ffv2\nLViw4PY79+/fn5qaGuGzpk6dOm7cuFZ3FhYWXmt70IS/xZ07q1Gj1IgRbR84Or16WfWVAeg0ebLa\nu1cdPqwCgQ4U0Js2bZoxY8aaNWuysrL8fn99fX1ZWVlNTU1RUVGEz4qLi4uLi2t1p9/vj/rMdefO\nqu2lG0BHlpqqOneW8zqhjoDOzs6uqKgoLi4+ffp0bW1tYmJiQUFBTk6O12vx6uE9MCFBRazqAHBT\naqpKSJDzR2N1BLRSyuv15ubm6lmrtawslZJiz9IAnCUzU6Wnqy++sHuOmzrAG1U2b1bx8XaPAsAJ\nhg1Tzz6rlJQXCV0d0GExHeAYAbQXSYkhaJT2F94DJX27AUjn8ShFg9aAgAYQrXBiENCaENAA2k5S\nYggapf2F98Dwf1gAoC1o0FpJ2g8BSCcpMQSN0v44Bw0gWjRorQhoAG0nKTEEjdL+aNAAokWD1oqA\nBtB2khJD0CjtjwYNIFq8UUUTAhpAtDjFoY9hCPlGA3AGSZVO0CjtzzBEfa8BOAANWh8CGkBUJIWG\noFHan2HwPm8A0aFB6yNpMwTgAJJCQ9Ao7Y9z0ACiRYPWhIAGEC0CWh8CGkBUJL1w5er8okEDiBYN\nWh8CGkBUJIWGoFHaHw0aQLRo0PoQ0ACiIik0BI3S/mjQAKJFg9ZH0guyABxAUqsTNEr7o0EDiBYN\nWhMCGkC0CGh9CGgAUZF0XtTV+UWDBhAtGrQ+BDSAqEgKDUGjtD8aNIBo0aD1IaABREVSaAgapf3R\noAFEiwatCQENIFoEtD6SrpgB4ACSWp2gUdofDRpAtGjQ+hDQAKJiGELSWbk8oGnQAEzweIRktNvz\ni4AGEC0xuSFlDkvQoAGYEBNDg9aCgAYQLTG5IWUOS9CgAZggJjekzGEJAhqACWLeP+H2/CKgAUSL\ngNbBMOR8owE4hphiJ2UOq4j5RgNwDDHFztX5xTloACYQ0JoQ0ACiRUDrQIMGYIKY3JAyhyUIaAAm\n0KA1IaABRIuA1oEGDcAEMbkhZQ6riPlGA3AMGrQONGgAJojJDSlzWEXMTgjAMcTkhqsDmgYNwAQC\nWhMCGkC0CGgdaNAATBCTG1LmsAQBDcAEGrQmBDSAaBHQOsTFqfh4u4cA4DQ+n90T3OTqgM7IUK++\navcQAJxmzx67J7jJ1QENAE5GQAOAUAQ0AAhFQAOAUAQ0AAhFQAOAUAQ0AAhFQAOAUAQ0AAhFQAOA\nUAQ0AAhlhEIhu2eIwoEDB5YuXdq1a9f7PvOzzz67ceOGYRgaptKssbExNjbW7inaX1NTk8fj4Ufm\nIC7+kcXFxWVkZLTlmXV1dcXFxSkpKZbMEXKpn/70p0eOHLF7ivZ36dKlWbNm2T2FJZYuXXrs2DG7\np2h/Fy9ezM/Pt3sKSyxZsuT48eN2T9H+Lly48NRTT9k9RSgUCnGKAwCEIqABQCgCGgCEIqABQCgC\nGgCEcm1Ax8TEeMT85cd2FBMTE+PSv4TLj8xx+JFZzWHXQbddQ0NDvEv/YqxbD82tx6Xce2huPS4l\n5tBcG9AA4HQiajwA4E4ENAAIRUADgFAENAAIRUADgFAENAAIRUADgFBuCOhPPvlk+PDhiYmJ8+fP\nb2hoiOpR4SIP//bbbw8cOLBLly4TJkwoKyuzZUJz2vJDKS8v9/l8mgd7cJEP7fz581OnTu3ateuo\nUaP+85//2DKhOZGP67XXXktPT09ISMjJySkvL7dlwgc0bdq0u05uc4DY/QupH1RjY2NKSsrmzZu/\n+OKL3NzcF154oe2PChd5+MrKSr/f/7e//e3q1au//OUvs7Ky7JozWm35oTQ1NY0ePdrj8egf70Hc\n99BGjhy5Zs2aysrKZcuW5eTk2DKkCZGP68SJE7GxsYcOHaqsrFyyZMmECRPsmtOcQ4cOLVy4UClV\nVlbW6iHbA8TxAX3o0KFBgwaFbxcXF/fv37/tjwoXefgdO3aMHTs2fLuhocEwjMuXL+se0ZS2/FDW\nrVuXn5/vuICOfGiffvrpwIEDg8FgKBQKBAIO+os/kY/r/Pnzfr+/tLS0rq5uxYoVjvuLPy+//PKS\nJUsSEhLuDGjbA8Sru7G3t1OnTg0ZMiR8Oysr68yZM8Fg8NYvOon8qHCRh580adKYMWPCt0tLS9PT\n07t3727PoFG67w/l1KlTGzZs2LNnT2FhoU0zmhT50I4ePTpgwICnn366uLh4yJAh69ats2/S6EQ+\nruTk5JdeemnUqFGGYfTo0cNZp26UUitWrFBK7dy5886HbA8QZ0RVBLW1tX6/P3y7a9euTU1N165d\na+OjwkUe3u/3P/TQQ6FQ6O233543b94rr7zilL/dGfm4gsFgQUHB2rVr2/KngaWJfGgXLlzYtWvX\niBEjdu/enZKSMnv2bJvGjFrk4yovL//1r39dUlJy/fr1hQsXLliwwKYx25/tAeL4Bp2YmFhfXx++\nXV9f7/F4bn9lKfKjwt13+JqamoKCgoqKip07d2ZnZ9sxoxmRj2vTpk0pKSnTp0+/dOmSTQOaF/nQ\nOnXqNH78+Kefflop9fLLL/t8vkuXLiUlJdkzazQiH9euXbsee+yxUaNGKaV+9atfdevW7erVq926\ndbNn1nZle4A4vkE//PDDty5gKC8vT09Pv/0/IJEfFS7y8A0NDZMnTx40aNCHH37ooHRW9zuuw4cP\nFxUVJSUlZWZmNjc3JyUllZaW2jRp1CIfWmpq6q3b4d847JRfphz5uJqbm4PBYPh2KBRqamoKueV3\nZNofIJrPebe78MusO3fuvHHjxsyZM1euXBm+f/v27efOnbvXo44Q+dC2bds2bNiwM7cJ/8OQL/Jx\nXbp06ezZs2fPnj1y5EhMTMzZs2cDgYCt80Yh8qHduHGjZ8+e27Ztu3r16nPPPTdu3Dhbh41C5OM6\nevRo9+7d33333StXrixbtmz8+PG2DmtS3759b3+RUEiAOD6gQ6HQRx99NHTo0J49e86fP//WP+Yu\nXbrs2rXrXo86RYRDe+6551rttdXV1fZO23aRf2Rh1dXVjruKI3S/Q3v//feHDh3q8/mmTp36+eef\n2zppdCIf11/+8pfMzEyfzzd9+vSKigpbJzWpVUALCRB+YT8ACOWYE7IA0NEQ0AAgFAENAEIR0AAg\nFAENAEIR0AAgFAENAEIR0AAgFAENAEIR0AAgFAENAEIR0AAgFAENAEIR0AAgFAENAEIR0AAgFAEN\nAEIR0AAgFAENAEIR0HAtr9fb1NS0b9++Tp062T0LYAYBDZcbOnToG2+8YfcUgBkENNxp8uTJzc3N\nGRkZlZWVq1evVkqdPHnykUceWbFiRVJS0tixY0tKSkaOHOn3+5955pnwp7z33nvf/OY3u3TpMnXq\n1MrKSlvHB5QioOFWBw4c8Hg8p06d6tKly607S0tLhw8f/t///jcQCMydO/ett97av3//unXrLl++\nXFNTk5eXt3r16nPnzmVkZHz/+9+3cXggzGv3AIA+ycnJ8+bNU0rl5uYqpdLS0tLS0vr27VtXV/fu\nu+9OmDBhxowZSqm1a9cmJSUFg8GYGBoM7ERAowPx+XzhG16v99Yrh16vVyl19uzZAwcOpKenh++M\njY29ePFinz597BgTuImABpRSqk+fPjNnzvzjH/+olAoGgxUVFb1797Z7KHR0/A8OblZfX9/GZ06b\nNm337t3FxcVXrlxZvXr1nDlzDMOwdDbgvghouFZ+fn6/fv2uX7/elienpKT86U9/Wrx4cUpKynvv\nvffnP//Z6vGA+zJCoZDdMwAA7oIGDQBCEdAAIBQBDQBCEdAAIBQBDQBCEdAAIBQBDQBCEdAAIBQB\nDQBCEdAAIBQBDQBCEdAAIBQBDQBCEdAAIBQBDQBCEdAAIBQBDQBCEdAAINT/A4QGkD2pILUiAAAA\nAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Optimal interval VWAP strategy in the square-root model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 2. (12 points)\n",
      "\n",
      "Consider equation (11) of Lecture 10 for the cost of liquidation using interval VWAPs:\n",
      "\n",
      "$$\n",
      "\\cC=\\sum_{j=1}^m\\,C_{jj}+\\sum_{i< j}^m\\,C_{ji}.\n",
      "$$\n",
      "\n",
      "(a) With $m=3$ (that is 3 interval VWAP slices), and assuming the square-root process where\n",
      "\n",
      "$$\n",
      "h(v)=\\frac  3 4 \\sigma\\,\\sqrt{\\frac{v}{V}} =\\frac  3 4 \\sigma\\,\\sqrt{\\frac{\\dot x}{V}}, \n",
      "\\quad \\mbox{ and } \\quad\n",
      "G(\\tau)=1/\\sqrt{\\tau},\n",
      "$$\n",
      "\n",
      "compute explicit expressions for the matrix elements $C_{ji}$ in terms of $x_1$, $x_2$, $t_1$, and $t_2$ (noting that $t_0=0$ and $t_3=T$).\n",
      "\n",
      "(b) Minimize the resulting expression for the expected cost of liquidation numerically to find optimal choices of $x_1$ and $x_2$ when $t_1=T/3$, and $t_2=2 \\,T/3$.\n",
      "\n",
      "(c) What is the percentage saving from trading with the 3-slice strategy relative to a VWAP?\n",
      "\n",
      "(d) Explain to your non-mathematical manager why this strategy has lower expected cost."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Solution\n",
      "\n",
      "### (a)\n",
      "\n",
      "For self-interaction term,\n",
      "\n",
      "$$C_{ii}=\\frac{x_i}{t_{i}-t_{i-1}}\\,h(v_i)\\,\\int_{t_{i-1}}^{t_i}\\,dt\\,\\int_{t_{i-1}}^{t}\\,\\frac{1}{\\sqrt{t-s}}\\,ds=\n",
      "\\frac{4}{3}\\,x_i\\,h(v_i)\\,\\sqrt{t_i-t_{i-1}} = \\sigma \\sqrt{\\frac{x_i^3}{V}} \\, \\forall i = 1 \\to 3$$\n",
      "\n",
      "\\begin{eqnarray*}\n",
      "C_{ji}\n",
      "&=& \\frac{x_j}{t_{j}-t_{j-1}}\\,h(v_i)\\,\\int_{t_{j-1}}^{t_j}\\,dt\\,\\int_{t_{i-1}}^{t_i}\\,\\frac{1}{\\sqrt{t-s}}\\,ds\\\\\n",
      "&=& \\frac{4}{3}\\,x_j\\,h(v_i)\\,\\frac{\\left(t_{j}-t_{i-1}\\right)^{3/2} -\\left(t_{j}-t_{i}\\right)^{3/2} -\\left(t_{j-1}-t_{i 1}\\right)^{3/2} +\\left(t_{j-1}-t_{i}\\right)^{3/2}}{t_{j}-t_{j-1}} \\\\\n",
      "&=& \\sigma\\,x_j\\,\\sqrt{\\frac{x_i}{V\\,(t_i - t_{i-1})}}\\,\\frac{\\left(t_{j}-t_{i-1}\\right)^{3/2}\n",
      "-\\left(t_{j}-t_{i}\\right)^{3/2}\n",
      "-\\left(t_{j-1}-t_{i-1}\\right)^{3/2}\n",
      "+\\left(t_{j-1}-t_{i}\\right)^{3/2}\n",
      "}{t_{j}-t_{j-1}}\n",
      "\\end{eqnarray*}\n",
      "\n",
      "Then we got,\n",
      "\n",
      "\\begin{eqnarray*}\n",
      "C_{21} &=& \\sigma x_2 \\sqrt{\\frac{x_1}{V}}(2\\sqrt{2}-2) \\\\\n",
      "C_{31} &=& \\sigma x_3 \\sqrt{\\frac{x_1}{V}}(3\\sqrt{3} - 4\\sqrt{2} + 1) \\\\\n",
      "C_{32} &=& x_3 \\sqrt{\\frac{x_2}{V}}(2\\sqrt{2} - 2)\n",
      "\\end{eqnarray*}"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### (b)\n",
      "\n",
      "The cost of liquidation,\n",
      "\n",
      "\\begin{eqnarray*}\n",
      "\\mathcal{C} &=& \\sum_{j=1}^m\\,C_{jj}+\\sum_{i< j}^m\\,C_{ji} = C_{11} + C_{22} + C_{33} + C_{21} + C_{31} + C_{32} \\\\\n",
      "&=& \\sigma \\sqrt{\\frac{x_1^3}{V}} + \\sigma \\sqrt{\\frac{x_2^3}{V}} + \\sigma \\sqrt{\\frac{x_3^3}{V}} + \\sigma x_2 \\sqrt{\\frac{x_1}{V}}(2\\sqrt{2}-2) + \\sigma x_3 \\sqrt{\\frac{x_1}{V}}(3\\sqrt{3} - 4\\sqrt{2} + 1) + x_3 \\sqrt{\\frac{x_2}{V}}(2\\sqrt{2} - 2)\n",
      "\\end{eqnarray*}\n",
      "\n",
      "where $x_1 + x_2 + x_3 = X$\n",
      "\n",
      "Without loss of generality, we let $\\sigma \\sqrt{\\frac{1}{V}} = 1$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "\n",
      "# Without loss of generality, we let the total number of shares equal to 1000\n",
      "X <- 1000\n",
      "\n",
      "liquidation_cost <- function(x) {    \n",
      "    x3 <- X - x[1] - x[2];\n",
      "    \n",
      "    cost <- x[1] * sqrt(x[1]) + x[2] * sqrt(x[1])*(2 * sqrt(2) - 2) + x3 * sqrt(x[1]) * (3 * sqrt(3) - 4 * sqrt(2) + 1) +\n",
      "        x[2] * sqrt(x[2]) + x3  * sqrt(x[2]) * (2 * sqrt(2) - 2) + x3 * sqrt(x3);\n",
      "\n",
      "    return(cost);\n",
      "}   \n",
      "\n",
      "##\n",
      "x_guess <- c(X/3, X/3)\n",
      "fitted_x <- optim(x_guess, liquidation_cost)$par\n",
      "cat(\"Given the total number shares is \", X, \". The x1 should be \", fitted_x[1], \" x2 should be \", fitted_x[2])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "Given the total number shares is  1000 . The x1 should be  647.0814  x2 should be  1.532216e-13"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "\n",
      "# (c)\n",
      "\n",
      "three_slices_strategy_cost <- liquidation_cost(fitted_x)\n",
      "VWAP_cost <- X * sqrt(X)\n",
      "cat(\"The percentage saving from trading with the 3-slice strategy relative to a VWAP is \", \n",
      "    (VWAP_cost- three_slices_strategy_cost)/VWAP_cost*100, \"%\")\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "The percentage saving from trading with the 3-slice strategy relative to a VWAP is  11.67179 %"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### (d)\n",
      "\n",
      "The reason is very straightforward. Because the market impact is transient. In this model, we are considering the market impact as transient which is true for the market. So, the cost will reduce as time goes by. But in the VWAP model, we assumed the market impact is permanent. The cost would not decrease as the time goes by. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}